Spectral formulation and WKB approximation for rare-event statistics in reaction 
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We develop a spectral formulation and a stationary WKB approximation for calculating the 
probabilities of rare events (large deviations from the mean) in systems of reacting particles with 
infinite-range interaction, describable by a master equation. We compare the stationary WKB 
approximation with a recent time- dependent semiclassical approximation developed, for the same 
class of problems, by Elgart and Kamenev. As a benchmark we use an exactly solvable problem of 
the binary annihilation reaction 2 A — > 0. 
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I. INTRODUCTION 

Since the pioneering works of Delbriick p|, Bartholo- 
may and McQuarrie et al. 0, 0, kinetics of react- 
ing systems with infinite-range interaction, containing a 
large but finite number of molecules or agents (such as 
bacteria, cells, animals or even humans) have attracted 
much attention 0, [(J . While the change in time of the 
average number of particles in such systems may be de- 
scribable by (continuum) rate equations, one often needs 
to know the probability of a non-typical behavior. This 
necessitates going beyond the rate equations, and a stan- 
dard way of achieving this goal is provided by the master 
equation of a gain-loss type which directly deals with 
n 2 ,...,n N (t)' the probability of simultaneously having 
n\ particles of the first type, na particles of the second 
type, . . ., and njv particles of the iVth type at time t 
. Though providing a complete description, the mas- 
ter equation is rarely solvable analytically, so various ap- 
proximations are in use 0, E| . Probably the most widely 
used approximation is the Fokker-Planck equation which 
usually suffices when one is not interested in extreme 
statistics, such as extinction time 0,0. Being interested 
in extreme statistics, we will exploit here a well-known 
mathematical formulation (see, e.g. 0) which, though 
exactly equivalent to the master equation, deals instead 
with a generating function G(x, t) which encodes all the 
probabilities P ni ,n 2 ,...,n N (t). The generating function is 
defined in the following way: 



G(x,i) = 



E 



"1=0, 
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whereas the probabilities P ni ,n 2 , 
differentiation: 



, (t) are recovered by 
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Til , 



1 d ni+ - +nN G(x,t) 
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This formulation transforms the master equation (an in- 
finite set of differential difference equations) into a single 
linear partial differential equation for G(x,t): 

<9G(x, t) 



dt 



LG{*,t), 



(3) 



where L is a real linear differential operator which in- 
volves derivatives with respect to the auxiliary variables 
x = (#1, X2, ■ ■ ■ , Xpj). The initial condition for this equa- 
tion is supplied by Eq. (yi with the time-dependent prob- 
abilities replaced by their (prescribed) values at t = 0. 
When the rate constants are time- independent, the op- 
erator L is independent of time. There is one universal 
boundary condition in this formulation. Indeed, as 



^ ] Pn 1 ,n 2 ,...,n N (t) — 1, 

n-i—0, riN—0 



one immediately gets 
G{ Xl = 1, 



,Xn 



l;t) = l 



(4) 



(5) 



Correspondingly, LG(x, t) must vanish at x = 

(X!,X 2 , ■ ■ -,X N ) = (1,1,..., 1). 

This paper will be limited to a single species, N = 1. 
In this case the generating function G{x, t) is 



G(x,t) = ^Tx n P n (t), 

n=0 

the probabilities P n (t) are 

1 d n G(x,t) 



Pn{t) 



dx n 



x=0 



the evolution equation for G is 

dG(x,t) 



dt 



LG(x,t). 



and the universal boundary condition (J^J is 

G(x = l,t) = 1. 
The initial condition is 



G(x,0) = G (x) = ^x n P n (Q) 



(6) 



(7) 



(8) 



(9) 



(10) 



n=0 



We will be interested in the important class of prob- 
lems where the operator L is of the second order [tj . The 
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crux of our approach is that, for any time- independent 
second-order linear differential operator L, one can ex- 
pand G(x, t) in a complete set of properly constructed or- 
thogonal spatial eigenfunctions of the operator L [IoL[ll| . 
The linear ordinary differential equation for these eigen- 
functions can be obtained, for any specific problem, by 
separation of variables. By a proper change of variables 
one can always eliminate the first derivative from this 
equation, and arrive at a spectral problem for a station- 
ary Schrodinger equation for a zero-energy particle in a 
potential V(x,fj.) which depends on a parameter \x com- 
ing from the separation of variables. This parameter, un- 
known a priori, represents the eigenvalue of this problem. 
This spectral formulation is exact, and it paves the way 
to a systematic computation of the probabilities P n {t), 
where n can be significantly different from the "typical" , 
or average number of particles. Using the probabilities, 
one can accurately estimate a host of quantities of inter- 
est, for example, the average extinction time and the life- 
time distribution. The present work is mainly concerned 
with one useful technique within the framework of the 
spectral formulation: the WKB approximation [12l ll.T ] 
which yields semiclassical eigenvalues and eigenfunctions. 
Using these, one can construct an approximate solution 
of the initial value problem for G(x, t) and, by virtue of 
Eq. 0, calculate P n (t) for n > 1. 

As Eq. ||3J| is readily interpretable as a (non- 
Hermitian) time-dependent Schrodinger equation with 
imaginary time, there were earlier quantum mechanical 
interpretations of rare event statistics in reacting sys- 
tems, such as the Doi-Peliti formalism of second quanti- 
zation 0, [TH fltl . Still another approach to this class 
of problems has been recently suggested by Elgart and 
Kamenev • Instead of dealing with the creation and an- 
nihilation operators, as is customary in the second quan- 
tization approach, Elgart and Kamenev reformulated the 
time-dependent problem in semiclassical terms, employ- 
ing two strong inequalities: n» 1 and n{t) 3> 1, where 
n(t) is the average number of particles in the system at 
time t. They showed that classical dynamics correspond- 
ing to the Hamiltonian of the problem provide a valu- 
able information about the rare-event statistics, and that 
this approach is greatly superior to the more customary 
Fokker-Planck description. Elgart and Kamenev start 
with the ansatz G(x, t) = exp[— S(x, t)] in Eq. (JSJ. Then, 
neglecting the d xx S term, they arrive at a Hamilton- 
Jacobi equation for S(x, t) which, for a time-independent 
L, is solvable [17]]. This procedure yields G(x,t) and, by 
virtue of Eq. J7J), the probabilities P n (t). Elgart and 
Kamenev illustrated this approach on several pedagogical 
examples which included various combinations of binary 
annihilation, branching, decay and creation of particles. 

As explained above, we suggest in this work a dif- 
ferent type of semiclassical approximation for this non- 
Hermitian quantum mechanics: a stationary WKB ap- 
proximation based on an exact spectral formulation. We 
will show that the two semiclassical approximations com- 
plement each other, each of them being advantageous in 



some region of the parameter space. We will demonstrate 
our approach by a simple example of binary annihilation 
reaction 2 A 0, where A > is the rate constant. In 
this case the master equation is 

j t P n {t) =^[(n + 2)(n + l)P„+ 2 (i) - n(n - l)P n (t)] , 
while the evolution equation for G(x, t) takes the form 



dG 
~dt 



X 



d 2 G 
dx 2 



(12) 



This example is instructive for two reasons. First, it is 
one of the examples used by Elgart and Kamenev 
to illustrate their time-dependent semiclassical approach. 
Second, and no less important, it is exactly solvable. Mc- 
Quarrie et al. p| used the exact solution to find the av- 
erage number of particles and the variance versus time, 
when starting from a fixed even number of particles. We 
significantly extend the analytical solution in Appendix 
and find the probabilities P n (t) for all n and t, and their 
various asymptotics. Using these findings, we also calcu- 
late the probability distribution of lifetimes of the parti- 
cles in this system and the average extinction time. These 
exact results provide a benchmark for our WKB theory. 
In principle, the strong inequality n 3> 1 is the only cri- 
terion required in this theory for all times. We observed, 
however, that in practice one also needs to require n > n, 
in order to avoid loss of accuracy resulting from summa- 
tion of many large terms of alternating sign. We will show 
that, in the region of n > n, the stationary WKB formal- 
ism is much more accurate than the time-dependent for- 
malism due to Elgart and Kamenev, while in the region 
of n < n the time-dependent formalism (which circum- 
vents the summation of large terms of alternating sign) 
is advantageous. 

Here is a layout of the rest of the paper. In Section 
II we apply separation of variables to Eq. I|12[l . arrive 
at a Sturm-Liouville eigenvalue problem and find ap- 
proximate solutions to this problem by using the WKB 
approximation. Then we solve, in Section III, an ini- 
tial value problem, compute the respective approximate 
probabilities P n (t) and compare them with the exact 
probabilities and their asymptotics, derived in Appendix. 
In Section IV we compare the predictions of the time- 
dependent semiclassical formulation with the exact 
results and establish the validity of the time-dependent 
and stationary semiclassical approximations. Section V 
presents a brief summary and discussion of our results. 



II. SPECTRAL FORMULATION AND WKB 
APPROXIMATION 

A. Boundary conditions and steady state solution 

To complete the formulation of the problem for Eq. 
1)12(1 . we need two boundary conditions. The first of 
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them, Eq. JJjJ, is universal. The second one, at x = — 1, 
readily follows from Eq. itself: G(x = — l,t) = C = 
G Q {x 



T) is determined by the ini- 



const, where G 
tial data fl8j . 

The limit of t — > oo corresponds to a steady state solu- 
tion of Eq. JEJl: G s (a;) = A + £?x. To obey the bound- 
ary conditions at x = ±1, we choose A = (1 + C)/2 
and B = (1 — C)/2. When C = 1 (an even number of 
particles at t = 0), the steady state solution G s (x) = 1 
corresponds to an empty system: P n = <5 Il0 , where <5y is 
the Kroenecker delta. When C = — 1 (an odd number of 
particles at t = 0), one obtains G s (x) — x. This corre- 
sponds to a single particle, P n — 5 n \, which lives forever 
as there are no particles it can react with. 



B. Separation of variables and eigenvalue problem 

Now let us consider the time-dependent part of the 
generating function: g[x, t) — G(x, t) — G s (x). As g(x, t) 
must vanish at x = ±1, we can look for solutions of the 
equation for g(x,t), 

dg 
dt 

in the separable form g(x, t) = exp(— jt) f(x). We obtain 



(13) 



<p»(x) 



1 — X 2 



0, 



(14) 



where /i 2 = 2j/X, and <£>(±1) = 0. Equation l|14|) can 
be interpreted as a stationary Schrodinger equation for a 
zero-energy particle (m = h = 1) in the singular potential 



V(x) = 



if 

2(l-x a ) ' 



N <i 

bl > 1. 



(15) 



see Fig. ^ with (a priori unknown) magnitude /i 2 which 
plays the role of eigenvalue. The problem is exactly solv- 
able in terms of Legendre polynomials, and the solution 
is presented in Appendix. In the next subsection we 
will proceed as if we were unaware of the exact solution, 
and find the spectrum of /i and the eigenfunctions in the 
WKB approximation. 



C. Stationary WKB approximation: wave 
functions and quantization 

The crucial assumption of our semiclassical theory is 
that the main contribution to the probabilities P n (t) with 
n 1 comes from the semiclassical region of spectrum of 
fi, where /i> 1 and the eigenfunctions have multiple ze- 
ros on the interval |x| < 1. We will verify this assumption 
a posteriori. Employing the str ong inequality /i» 1, we 
find by a standard calculation 0, ^ two independent 
(even and odd) WKB solutions of Eq. l (T4*)l : 



Revert yE) 
Vodd{x) 



(1 
(1 



,2)1/4 
r 2 W4 



cos(/iarcsina;) 
sin(/i arcsinx) . 



(16) 
(17) 



£0 



E=0 



-1 





X 



FIG. 1: (Color online) Shown by the blue solid lines is the 
singular potential V(x), given by Eq. I|15|l . The vertical scale 
is arbitrary. 



To quantize the eigenvalue /i, we must use the boundary 
conditions at x = ±1. The WKB solutions (|16|) and 
(|17f) are invalid, however, at and near the singular points 
x = ±1. To determine the solutions there we first solve 
Eq. (|14fl in a small vicinity of each of these points. Here 
it suffices to consider the point x = 1. Let us introduce a 
new coordinate £ = 1 — x -C 1 and neglect the subleading 
terms in the expansion of the potential V(x) near x — 1. 
Equation (|rl|) becomes 



The two independent solutions are 

pi(o=e 1/2 -wo 1/a ] 

and 

M0 = e /2 YM20 1/2 }, 



(18) 



(19) 



(20) 



where J\ and Y\ are the Bessel functions of the first and 
second kind, respectively. Only <^i(£) obeys the required 
boundary condition tp(£ = 0) = 0, so ^a(C) m ust be 
discarded. Now we can match </?i(£) with each of the 
WKB solutions ||T^| and JTJJ. Indeed, when fi > 1, 
the solution <^i(£) remains valid at /i^ 1 ^ 2 ^> 1, as long as 
£ <C 1. It has, therefore, a common region of validity with 
the WKB solutions. We use the asymptotic expansion 

M 



Ji(z) ~ W— cos — 

V 7TZ \ 4 



at z > 1 



(21) 



and obtain, up to a constant factor, 

^(0 ~ e^(^-l 



(1 - x) 1 /* sin \^yj2(l-x) - j 



(22) 
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Now we expand the even WKB solution (|16fl at 1— x <C 1: 

<p ev (x) ~ (1 - x) 1 ' 4 sin [ Mv /2(1 - ^ - |] . (23) 

Matching the asymptotes (1221) and i|23[l . we obtain the 
discrete spectrum eigenvalues corresponding to the even 
eigenfunctions: 



/i = 2k 



1 

2 ' 



(24) 



where fc 3> 1 is an integer. In a similar way we obtain 
the discrete spectrum for the odd eigenfunctions: 



(25) 



The eigenvalues (|24|l and l|25|) coincide, in the leading 
and subleading orders in k 3> 1, with the exact eigen- 
values (4fc 2 — 2A:) 1 / 2 and (4k 2 + 2k) 1 / 2 , respectively, see 
Appendix. Furthermore, the corresponding WKB eigen- 
functions provide an accurate approximation, see Fig. [21 
to the exact eigenfunctions [which are orthogonal, with 
respect to the inner product with the weight function 
w(x) = (1 — a; 2 ) -1 , on the interval < 1, and form a 
complete set] . We normalize the approximate eigenfunc- 
tions ui(x) by demanding 



u 2 (x)w(x) dx 



The normalized even WKB eigenfunctions are 



(26) 



-2\l/4 



x cos [fik arcsin x] 

for < \x\ < 1 -£, 

4^(x) =-V2/i fe (l-N) 



(27) 



X Ji 



^2(1- M) 

for 1 — e < \x 



< 1. 



where 1/fc 2 < e < 1, ^fe = 2k - 1/2, and fc > 1. For 
l//i 2 -C 1 — |x| -C 1 the function u^i (%) coincides, in the 



leading order, with u2 (&). It is sufficient to use only the 
v,2 k {x) asymptotes in the normalization integral (|26J) . 



III. INITIAL VALUE PROBLEM AND 
CALCULATION OF THE PROBABILITIES 

Putting everything together, we can write the WKB 
solution of the initial value problem for Eq. (|12|l as 



G(x,t)~G s {x) 



l>0 



ai exp 



ui(x), (28) 



where each constant a; is equal to the inner product 
[with the weight function w(x) — (1 — a; 2 ) -1 ] of the re- 
spective normalized eigenfunction ui (x) and the function 
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FIG. 2: (Color online) The normalized WKB eigenfunctions 
127H (the red dashed lines) and the normalized exact eigen- 
functions (the blue solid lines) for k = 10 (a) and k = 2 (b). 
As one can see, a good agreement is observed even for k — 2, 
while for k = 10 the agreement is excellent. The inset in (b) 
shows a close vicinity of x — 1. One can see that the function 
u 20 ( x ) (the red dashed line) deviates from the exact solution 
(indicated by the crosses) in the vicinity of x — 1, whereas 
the function itjo' ( x ) (the solid line) is in excellent agreement 
with the exact solution there. 



Gq(x) — G s (x). One can see that the populations of the 
eigenstates with I > of this non-Hcrmitian "quantum 
mechanics" are decaying exponentially in time. 

Assume, for concreteness, that the initial number of 
particles is fixed and equal to uq = 2ko 3> 1, where ko 
is integer. In this case [see Eq. <|10[l ] Go(x) = x 2k °, and 
one only needs the even eigenfunctions l|27|) . Now we can 
use the orthonormality relation l|26|l and compute the 
coefficients a k ■ After a lengthy algebra we obtain 



a h . 



2 fc(2fc-l) 

. e 2k o 



(29) 



where we have assumed k < y/k~o <C fco, as justified below. 
In the leading order in 1/fc < 1 Eq. Q29JI coincides with 
the corresponding asymptotics (|A4|I of the exact result. 
Now Eq. J2HI becomes 



G(x,t) ~1 + 



OO 

E 

k=l 



-fc(2fc-l) 



~u 2 k(x) ■ 



(30) 
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where 



N(t) = 



2fc 
1 + 2k \t 



is the average number of particles at time t according to 
the mean-field theory. Though the sum in Eq. (|30[l for- 
mally runs to infinity, the dominant contribution comes 
from terms with k < \fk~o <C fco i while the rest of terms 
give only exponentially small corrections. 

To recover P n (t), for 1 <C n < V^Oi we use Eq. (JJJ: 



Pn(t) 



1 



OO 

E 

fc=i 



d n U 2 k(x) 

g W(t) ZKV ' 



9: 



(31) 











-1 









-1 






1 
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3 > 





As our WKB approximation assumes fik = 2k — 1/2 ^> 1, 
the first few terms of the sum in Eq. 1)31(1 may seem 
inaccurate. It turns out, however, that the sum actually 
starts from k = n/2, see below. Therefore, at 1 <C n < 
V^fco all the terms of the sum in Eq. Il-illl are accurate. 

Equation l|31|) is the central result of the stationary 
WKB approximation for the binary annihilation prob- 
lem. To compute the n-th derivative of the WKB eigen- 
functions u 2 k{x), entering Eq. I|31l) . we can analytically 
continue u 2 k(x) into the complex plane. By virtue of the 
Cauchy theorem 



I = 



1 d n u 2 k(x) 



0: 



1 

2nd 



u 2k {z)dz 



(32) 



where the integration is performed over a closed contour 
C in the complex z-plane around the pole x = inside the 
region of analyticity of u 2 k(x). We can put here u 2 k{z) — 

(z), since the main contribution to the integral, as 



U 

shown below, comes from the region far from the points 
x = ±1, in the vicinity of which ttfc (z) is inaccurate. We 
obtain 



(~l) fc 
V2tt 3 / 2 « 



Im 



g{z)e {n ^ f{z U2 



G 



(33) 



and 



where f(z) = iAarcsm(z) — ln(z), g(z) 
A = /ifc/(n — 1). As n 3> 1, the integral in (|33|l can 
be evaluated using the saddle point approximation [l3| . 
This is done by deforming the contour C, so that it 
passes through the saddle point z* = x* + ij/*, where 
u(x,y) = Re[/(z)] obtains its maximum, and conse- 
quently v(x,y) = Im[/(z)] is constant. By choosing 
the contour at the saddle point to be parallel to the 
direction of the steepest descent of u(x,y), we can re- 
place the integration over the complex plane by inte- 
gration over the real axis, having to multiply the result 
by a constant phase: the value of v(x, y) at the saddle 
point. The saddle point can be found from the equation 
f'(z) = 0. For k < n/2 (that is, A < 1), the saddle point 



FIG. 3: A possible contour C, see Eq. . The direction 

(2) 

of the contour in the vicinity of the saddle point z — zi 
(marked by the circle) is chosen to be along —X7u(x, y) at 
Zi , implying in this case a — 0. As n > 1, the contribution 
of the rest of the contour (as long as it encircles z — 0) is 
exponentially small. 



Z * 2) = 

becomes 113 



-1/2 



In each of these cases Eq. (|33f) 



I = 



(-1)* 



Im 



2ng(z* 



= (n-l)/(z„). 



V(n-l)|/"(z,)| 



(34) 



where a is the angle of the contour with respect to the 
positive real axis at the saddle point, where the contour 
is chosen to be parallel to — S7u(x,y). 

This procedure yields markedly different results in the 
cases of A < 1 (k < n/2) and A > 1 (k > n/2). For 

A < 1, upon substituting z^ and deforming the contour 
so that a = 7r/2 in its vicinity, we realize that the result 
inside the brackets in Eq. (|34|l is real which yields 1 = 0. 
Therefore, the saddle-point asymptote Eq. (|34|l predicts 
that the n-th derivative of the eigcnfunctions u 2 k (x) van- 
ishes for k < n/2, so the sum in Eq. H31|) starts from 
k = n/2. This could be expected, as the same kind of 
behavior is exhibited by the exact eigenfunctions q 2 k (x) , 
which are polynomials of order 2k, see Appendix. 
After some algebra, Eq. (|34|l yields, for k > n/2 



( _ 1) fc-„/2 e l/4 A .3/2 2 „ (fc + |) 



(35) 



where we have substituted zi ; for the saddle point. A 
possible contour C is shown in Fig. |3| We can see now 
that, as the saddle point z* lies on the imaginary axis, 
the contour does not have to come close to x = ±1, thus 
justifying the use of u$(z) for u 2 k(z). The saddle point 
approximation is only valid when \ f" {z*)\/\f" {z*)\ z / 2 <C 
[1 — A 2 ]~ 1 / 2 , whereas for ^Jn 01 . For k > n/2 this requirement is equivalent to 

k — n/2 ^> 1. We will see shortly, however, that the 



lies on the real axis: 



k > n/2 (that is, A > 1) it lies on the imaginary axis 



results remain quite accurate even for k — n/2 — 0(1) 
Now we can rewrite Eq. (|31fl as 



Pn(t) 



fe+#-l 



7r n 

R «•(*) 



(36) 



We immediately notice that this formula is very similar to 
Eq. (|A13|) . Moreover, the two expressions coincide if we 
rewrite the factor (k — n/2) k ~ n / 2+1 / 2 in the denominator 
of Eq. (|XT3)) as (k - n/2 + 1/4 - i/4)fe-"/2+i/2 ) assume 
fc — n/2 ^ 1 and use the asymptote (1 + £/u) u ~ at 
« ^> 1. As Eq. I|A13() gives an accurate approximation 
to the exact probabilities at n ~ TV (see Appendix), the 
same is true for the stationary WKB result l|A13[l . For 
n ^> N it suffices to take into account only the first term, 
k = n/2 in the sum of Eq. 136fl which yields 



Pn{t) 



2 e<i 

g 2JV(«) 



(37) 



This asymptote coincides, up to a factor y/2/ne 1 ^ 4 ~ 
0.976, with Eq. I|A12(I that gives an accurate approxima- 
tion to the exact probabilities in this limit, see Appendix. 
Therefore, in the region of n > N the stationary WKB 
theory is accurate, as can be seen from Figures 14161 Im- 
portantly, the asymptote l|37|) does not demand TV ^> 1 
and therefore remains valid at long times, Xt £ 1, when 
the average number of particles is already small, see Fig. 
□ 



IV. TIME-DEPENDENT SEMICLASSICAL 
SOLUTION VERSUS EXACT SOLUTION 

The time-dependent semiclassical approximation, sug- 
gested by Elgart and Kamenev 0, differs from the sta- 
tionary WKB approximation in that it deals semiclassi- 
cally with the original non-Hermitian Schrodinger equa- 
tion JSJl (a partial differential equation) , rather than with 
the set of ordinary differential equations obtained by the 
separation of variables in Eq. (JSJl. Using the ansatz 
G(x,t) = exp[— S(x, t)] and neglecting the d xx S term, 
one arrives at a Hamilton-Jacobi equation. For the bi- 
nary annihilation problem this equation is 



dS X , 2 , fdS 



(38) 



Elgart and Kamenev found an exact solution to this 
equation: 

S(x, t) = - N(t) arccos 2 x . 
The respective generating function 



G(x, t) = exp 



- N(t) arccos 2 x 



(39) 




FIG. 4: (Color online) The decimal logarithm of P„(t) as a 
function of n/N for Xt = 0.02 and no = 10 3 . Shown are 
the exact probabilities (IA8II (the solid line), the stationary 
WKB probabilities l-'i U (the empty circles) and the time- 
dependent WKB probabilities 1 1 U (the blue squares). In- 
set shows the logarithm of the ratio of the stationary WKB 
probabilities and the exact ones (the red solid line), and the 
logarithm of the ratio of the time-dependent WKB probabili- 
ties and the exact ones (the blue dashed line), versus n/N. As 
n/N grows, the time-dependent WKB solution deteriorates, 
whereas when n goes down below N, the stationary WKB 
solution deteriorates. 
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FIG. 5: Same as in Fig. H but for Xt = 0.1. 



obeys the boundary condition Q exactly, and the initial 
condition G(x, t = 0) = x n ° with a high accuracy as long 
as no 3> 1. The probabilities P n (t) can be calculated 
from the equation 



Pn(t) 



1 d n G(x,t) 



dx r - 



x=0 



±J^G( Z ,t), (40) 



where the integration is performed over a closed contour 
including z = in the complex z plane, inside the region 
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FIG. 6: Same as in Fig. but for At = 0.5. One can clearly 
see that the time-dependent WKB result is already inaccurate 
at n > N , while the stationary WKB result keeps its high 
accuracy there. 
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FIG. 7: Long time asymptotics of P n (t). Shown is the decimal 
logarithm of P„(t) as a function of At for no = 10 3 and n = 
2, 4, and 6. The circles denote the stationary WKB solution 
[Eq. 1371 1. the solid lines denote the exact solution [Eq. (IA8H . 
Notice excellent agreement even for small n. 



of analyticity of G(z,t). For n > 1 and N(t) > 1, the 
integral can be evaluated by the saddle point approxima- 
tion, and the result is 



Pn(t) 



1-x*. 



2tt ra - N{t)x 2 s 



; -JV(t)(| arccos 2 x s + jffe lni a ) 

(41) 



where x s = x s (n/N) is the root of the saddle-point equa- 
tion x s (l — ij)" 1 / 2 arccos:r s = n/N j7(- For n <C N(t) 
one obtains x s ~ 2ti/(ttN), so 

, 7r7V(t) 7T 2 7V(t) , s 

lnP n (i)~nln-^ ^ + n (42) 



(in the corresponding asymptote of Ref. the term 

A 

3[n-N(t)} 2 



is missing.) 

For n ~ 2V(f), x a ~ 1 - (3/2) (1 - n/iV), so 



lnP„(f) l- 

Finally, for n > iV(t), x s ~ (1/2) e n /^, so 
lnP„(t) ~ rain 2 - 



2iV(i) 



(43) 



(44) 





FIG. 8: The dashed regions schematically show the validity of 
each of the two semiclassical theories on the parameter plane 
(N, n). Figure a shows the validity domain n > N and n> 1 
of the stationary WKB theory. An additional condition for 
its validity is n < \fk~o <C ko (not shown). Figure b shows 
the validity domain 1 <C n < N of the time-dependent WKB 
theory. 

What arc the applicability conditions of the time- 
dependent semiclassical approximation? The above cal- 
culations required n = 2k 3> 1, n ^> 1 and N ^> 1. To 
find out whether there is an additional condition, let us 
consider the n 3> N(t) asymptote of the exact result [Eq. 
(EH)]: 



lnP n (i) ~ nln.2- 



2N(t) 2N(t) 



(45) 



A comparison of Eqs. H44|) and (|45H shows that the time- 
dependent WKB probability P n (t) lacks a large term 
n/[2N(t)] in the exponent, and therefore it greatly un- 
derestimates rare events with n 3> N(t). This effect can 
been seen in Figs. 14161 On the other hand, in the region 
n < N(t), the time-dependent WKB theory yields a good 
approximation to the exact result, as it circumvents the 
summation of large terms of alternating sign. 

Therefore, based on the analytical and numerical com- 
parisons, we conclude that the time-dependent semiclas- 
sical approximation is accurate at 1 <C n < N. This 
obviously implies iV > 1, that is not too long times: 
At <C 1. The stationary WKB approximation is accurate 
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for n> N for any N, that is for all times. These results 
are illustrated in Figs. I4l6l which show P n (t) versus n/N: 
the exact result l|A8ll and the predictions of each of the 
two approximations, Eq. (|31|l and Eg. I|41l) . The validity 
domains of each of the two semiclassical approximations 
in the parameter plane (N, n) are shown in Fig. |H1 



gral of motion, the classical motion is non-separable and, 
in general, chaotic |2(j. Therefore, for highly-excited 
states, where classical mechanics and stationary WKB 
approximation are relevant, the spectral formulation may 
bring about an extension of quantum chaos to these non- 
Hermitian "quantum" systems. 



V. SUMMARY AND DISCUSSION 
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We developed a spectral formulation and a stationary 
WKB approximation for calculating the probabilities of 
rare events in systems of reacting particles with infinite- 
range interaction which are describable by a master equa- 
tion. We extended the exact analytical solution of the 
binary annihilation problem 2A — > and used is as a 
benchmark for testing the stationary WKB approxima- 
tion and a recent time-dependent WKB approximation 
due to Elgart and Kamenev Q. In theory the station- 
ary WKB approximation is always more accurate than 
a time-dependent WKB approximation. In practice this 
advantage is indeed realized in the regimes where the 
superposition of the different "quantum states" of the 
system is dominated by a small number of terms. On 
the contrary, when many "quantum states" are involved, 
virtually any approximation to the "wave functions" of 
individual states may alter the precise destructive inter- 
ference between the different states and cause large er- 
rors. In such cases the time-dependent WKB approx- 
imation 0, which effectively sums over the "quantum 
states" without dealing with them explicitly, can be ad- 
vantageous. 

WKB approximation alone is insufficient for calcu- 
lating the life-time probability distribution of systems 
which exhibit extinction. This quantity is encoded in 
Po(t) = G(x = 0, t) which involves a sum over all "quan- 
tum states" , including the lowest ones (see the final part 
of Appendix) . Still, the spectral formulation can be very 
useful here: it clearly identifies the lowest states and 
provides a proper framework for calculating their eigen- 
values and eigenfunctions: for example, by a variational 
method. 

This work dealt with the case when the operator L 
is of the second order, and the standard machinery of 
Sturm-Liouville theory is therefore available. The spec- 
tral formulation itself, however, is equally applicable to 
higher-order operators. Furthermore, it is well known 
that "WKB analysis is not sensitive to the order of a dif- 
ferential equation" (Ref. p. 496) which paves the 
way to generalizations of the theory to more complicated 
reaction kinetics. 

Finally, we have restricted ourselves in this paper to a 
single species. The generating function formalism, how- 
ever, is applicable to any number of species [see Eqs. 0- 
Already for two species some qualitative changes are 
possible. Indeed, for two species, the underlying clas- 
sical phase space, described by the Hamiltonian of the 
problem, is four-dimensional. If energy is the only inte- 
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Appendix 

Here we briefly review and extend the exact solution of 
the binary annihilation problem, obtained by McQuarrie 
et al. 4]. Let the initial state correspond to a fixed 
and even number of particles, so one only needs the even 
eigenfunctions of Eq. (|14|l . with the boundary conditions 
cf>(x — ±1) = 0. The exact even eigenfunctions are [2l| 

<f>2k{x) = T>2k{x) - T>2k~2{x) , k = 1, 2, . . . , 

where Vi(x) is the Legendre polynomial of order I. The 
corresponding exact eigenvalues are /i = (4fc 2 — 2k) 1 / 2 . 
The normalized even eigenfunctions, see Eq. 1)26(1 . are 



l2k{x) = 



k(2k- 1) 
4k- 1 



[V 2k (x) - V 2k -2(x)} 



The exact solution of an initial value problem for G{x, t) 
can be written as 



G(x,t) 



E 

k=l 



A kq2k {x)e- k ^ k -^ 



(Al) 



where the coefficients A k are determined by the initial 
data Go(x). For the initial data we are interested in no = 
2k , where k = 0, 1, 2, therefore, Gq(x) — x 2k " . Using 
the orthogonality relations of the Legendre polynomials, 
we obtain, after some algebra, 



A, 



4k -1 



r(i + k )r (i + fe ) 



k(2k - i) r (| + fc + k) r(k -k + i) 



(A2) 



where T(- ■ •) is the gamma function. Owing to the pres- 
ence of factor T(ko — k + 1) in the denominator, A k van- 
ishes for k > k . Therefore, the sum in Eq. i|Al(l is finite 
in this case, and it ends at k = fco- The average num- 
ber of particles at time i, n(t), can be found from the 
following relation: 



n(t) 



71 = 



nP n {t) 



dG{x,t) 



Ox 
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Using Eq. (|A1|) . we obtain 



i(t) = ^/k(2k- l)(4k- 1) A 



fee 



-k(2fc-l)A* 



(A3) 



fc=i 



Equations (|Al |l -(jA3 p coincide, up to notation, with the 
results of McQuarrie et. al. 4| [they also calculated the 
second moment n 2 (t)\. We now extend the exact theory 
in three directions. First, we obtain some useful approxi- 
mations for a large number of particles. Second, we calcu- 
late the probabilities P n (t) and their approximate asymp- 
totics in different regimes. We use these approxima- 
tions while comparing the exact solution (i) with our sta- 
tionary WKB results, and (ii) with the time-dependent 
semiclassical approximation of Elgart and Kamenev Q. 
Third, we find the probability distribution of lifetimes of 
the particles in this system, and the average extinction 
time. 

When k < \fk~Q (see below), we can reduce Eq. (IA2|I 

to 



4k- 1 



y fc(2fc — i) 

which yields the following approximation for G(x, t) 



(A4) 



G(x,t) 



E 



k=l 



I 4k- 1 
k(2k - 1 



k(2k-l) 

«(*) . (A5) 



Here 



N(t) 



2k 



1 + 2k Xt 



which is the mean-field result for the average number 
of particles. Equation (|A5(I is valid for all times, as the 
dominant contribution to the sum comes from terms with 
k \k~o ^ ^Oi while the rest of terms give only exponen- 
tially small corrections. The average number of particles 
(|A3II can be approximated as 



n(t) 



k 

E 

k=l 



(4k - 1) e 



(A6) 



For short times, Xt <C 1, N(t) 1, so the summation in 
Eq. I|A6() can be replaced by integration. Moving the up- 
per limit to infinity (which only causes an exponentially 
small error), we obtain 



J k(2k-l) 

(4k - 1) e "W dk 



N(t), 



(A7) 



the mean-field result. In the long-time limit Xt 1, 
N(t) ~ (Xt)- 1 < 1, the term k = 1 in Eq. (JXBJ is 
dominant, and we obtain n(t) ~ 3e~ At . 

Now we employ Eq. JJJ of the main part of the paper 
to calculate the probabilities P n (t). After some algebra 
we obtain the exact result: 



P n (t) = 6 0n + 



2 n- 



J2C e (k,n)e^ 2k ' 1 



(A8) 



where n is assumed to be even, r = max(l,n/2), and 
"(-i) fe ~ f (4fc-l)r[£:-i + ! 



C e (k,n) 



A, 



0F(fc-3)! 



(A9) 



At 7i > the sum in Eq. (| Afi|> starts from k = n/2. 
This is because the eigenf unctions qik (x) are polynomials 
of order 2k, so the n-th derivative of q2k(x) vanish at 
n > 2k. 

Going to the limit of n 1 and k < \fk~o : and using 
the approximation ljA4fl for A^, we can rewrite C e (k,n) 
as 



C e (k,n) 



(-l) k -^2 5 / 2 k(k + |) 



fc+S-i 



fc(2fe— 1) 



(*-$) 



Therefore, at 1 <C n < V^Oj -fn(*) becomes 



(A10) 



2 n+l ^ (-l)fc-ffc(fc+|) fc+ ^ 1 -K2*-i) 

(All) 

This expression can be simplified drastically for n 3> 
N(t). Here the sum can be accurately approximated by 
its first term k = n/2, and we obtain 



Pn(t) 



-n(n-l) 

e 2Jv w 



(A12) 



This asymptote coincides, in the limit of 1 <C n % \fka, 
with the first term of the sum in the exact result 




k—r 



FIG. 9: The exact probabilities (IA8I (the solid line), the ap- 
proximate probabilities 1A111 (the circles), and their asymp- 
totics (|A 1 3|> and l|A12fl (the squares), for no = 10 4 and 
N(t = 0.01) ~ 99.5. The agreement is good for n > N(t). 



In the cases of n <C N(t) and n ~ N(t) the leading 
contribution to the sum in Eq. (fATTTl comes from terms 
for which k — n/2 3> 1 which makes it possible to use the 
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Stirling formula for the factor (k — n/2)\ and arrive at 



Pu{t) 



•>n+i 



k 

E 

*;>» 



(-l) fc -7fc(/c + f) 



k+i 



-fc(2fc-l) 

■ e 



(A13) 



Now let us go back to Eq. (jAlip . According to our 
analysis, it should be accurate for 1 -C n < Vfa). A nu- 
merical comparison with the exact result (|A8(I [see Fig. 
[5] shows, however, that Eq. IjAllfl is accurate only at 
N < n < \Ao- At n < A the agreement rapidly dete- 
riorates. The disagreement stems from the fact that the 
sum in Eq. (|A8|) consists of terms of alternating sign. In 
the region of n < A, P n (t) is much smaller than each of 
the relevant terms of the sum, while the magnitudes of 
the successive terms are close to each other. One can say 
that there is strong destructive interference of "quantum 
states" of the system. In this situation, virtually any 
approximation made in calculating the individual terms 
of the sum may alter the precise balance between the 
terms and cause large errors in the region of n < A. The 
same problem appears in our stationary WKB theory, 
see Section III, which makes the time-dependent WKB 
approximation advantageous in this case. 




FIG. 10: (Color online) The lifetime probability distribution 
p(t) normalized to A, as described by Eq. 1A15I (the blue 
solid line), and its long time asymptote p(t)/\ = (3/2) e~ xt 
(the red dashed line). 



Now we proceed to calculating the probability distri- 
bution of lifetimes of the particles, and the average ex- 
tinction time. The quantity Po(t) is the probability of 
extinction at time t. Therefore, the lifetime probabil- 
ity distribution is p(t) — dPo(t)/dt. On the other hand, 
P (t) = G(x = 0,t). Therefore, using Eq. we 
obtain the exact result 



p(t) = \J2k{l~2k)A k q 2k (0)e 



-k(2k-l)\t 



(A14) 



k=l 



Both p(t) and all its derivatives with respect to t vanish 
at t = 0, so pit) is exponentially small at Xt <C 1. When 
k 3> 1 and Xt > 1/ y/ko, Eq. (| A14|) can be approximated 



p(t) ~ A k(2k - 1) [P 2k -2(0) - V2k(0)} e 



k=l 



k(2k-l)\t 

(A15) 

which is independent of ko . This universal distribution is 
shown in Fig. The long-time tail of the distribution, 
At ^ 1, is described by the first term of the sum in Eq. 
(lATsll . which yields p(t) ~ (3A/2) e" At . 
The average extinction time is 

/•OO 

tp(t)dt = [1 - P {t)] dt = 

h] Jo 

[1 - G[x = 0,t)} dt, (A16) 



which yields 



-, kg 

E 



A k q2k(0) 



(A17) 



A ^ k(2k - 1) 
When fco 1, we obtain a fco-independent asymptotics 
1 ^ V2k(0) - V 2 k-2(0) _ 1.38629.. 



A ^ k(2k - 1) 



A 



(A18) 



which also follows from Eq. (|A15(1 . Note that the first 
term in the series of Eq. i|A18(l already gives a fair accu- 
racy: fi = 1.5/A. 
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